#!usr/bin/perl -w
use strict;

=head1 NAME

searchforTIR_TcMariner - find TcMariner TEs based on TIR and TSD extraction

=head1 SYNOPSIS

 perl searchforTIR_TcMariner.pl -db zoe_input_5_Nc_10.fa -t 5_vs_4_ncbi.blastn.filterd.tab

=head1 DESCRIPTION

 Provide a database of the query sequences and the BLAST result of a cannonical TE against the DB

=head1 AUTHORS

Zoe Zhou
Jason Stajich

=cut

use Bio::DB::Fasta;
use Bio::SeqIO;
use Bio::Perl;
use Bio::Tools::Run::WrapperBase;
use Getopt::Long;

my $Prefix = 'TcMariner';
my $Min_PID = 90; # require a minimum percent identity
my $TIR_max_search = 100;
my $TSD_max_search = 10;
my $TIR_min_length = 9;
my $TSD_f = 'TA';
my $TSD_r = reverse($TSD_f);
my $TSD_length = length($TSD_f);

my $wrapper = Bio::Tools::Run::WrapperBase->new;
my $bl2seq = $wrapper->executable('bl2seq'); # find bl2seq if it is in your path
my $db;
my $table;
my $noblast = 0;
my $outdir = 'Full_TcMariner';
GetOptions('d|db:s'    => \$db, # FASTA database of sequences from the BLAST
	   't|table:s' => \$table, # BLAST report in -m8 or -m9 format (or -m 8C from the FASTA program)
	   'o|outdir:s' => \$outdir, # place to store the output files
	   'bl|blast:s' => \$bl2seq, # if you know the location of the bl2seq program, and it is NOT in your PATH
	   'n|noblast!'   => \$noblast,
);

$table ||= shift @ARGV;
mkdir($outdir) unless -d $outdir;
open(my $fh => "<$table")  || die "cannot open file $table: $!\n";

#my $tmp = $db;
# strip the ending which can be several different things
#$tmp =~ s/\.(fasta|fa|seq|fas|nt)$//;
#my $tmp1 = $tmp;
#my $output = "$outdir/$tmp.full.fasta";

my $dbh = Bio::DB::Fasta->new($db);
my %data;
my $ct = 0;
while (<$fh>) {
  # this is parsing the BLAST file which has these columns

  # Fields: query id, subject ids, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score

  next if /^#/;
  chomp; # drop the ending "\n" since we are using '\t' as the split option
  my ($qid,$subjectid, $percentID, $aln_length, $mismatches, $gap_open,
      $qstart,$qend,$sstart,$send, $evaluem,$bits) = split(/\t/,$_);

  # I want to capture all the lines from the BLAST report
  # push @{$data{$qid}->{$subjectid}}, [$percentID,$aln_length,$qstart,$qend,$sstart,$send];

  ($qstart,$qend) = sort { $a <=> $b } ($qstart,$qend); # insure start < end
  my $seqlen = $dbh->length($qid);
  my ($left_s, $left_e) = ($qstart - $TSD_max_search, $qstart + $TIR_max_search);
  my ($right_s, $right_e) = ($qend - $TIR_max_search,$qend + $TSD_max_search);

  # deal with boundaries - if we go before start of seq or if go into space of the other sequence
  $left_s = 1 if $left_s < 1;
  $left_e = $right_s -1 if $left_e >= $right_s;

  $right_e = $seqlen if $right_e > $seqlen;
  $right_s = $left_e + 1 if $right_s <= $left_e;

  my $leftend = $dbh->seq($qid,$left_s => $left_e);
  my $rightend = $dbh->seq($qid,$right_e => $right_s); # ask for this is revcom coordinates by reversing these

  #######################################matching section beginning###############################################
  #this finds the longest common substring in the left and right end
  # can we replace this with an alignment based approach?
  my $longest_com_sub;
  my $length_TIR =0;
  my ($l_pos,$r_pos); # position of the TIR

  #warn "left string is  $leftend\n";
  #warn "right string is $rightend\n";

  if (  defined $bl2seq && ! $noblast) {
    my ($f1,$tf1) = $wrapper->io->tempfile;
    my ($f2,$tf2) = $wrapper->io->tempfile;
    Bio::SeqIO->new(-fh => $f1,-format=>'fasta')->write_seq(Bio::Seq->new(-seq => $leftend,-id => 'left'));
    Bio::SeqIO->new(-fh => $f2,-format=>'fasta')->write_seq(Bio::Seq->new(-seq => $rightend,-id => 'right'));
    close($f1);
    close($f2);
    open(my $run => "$bl2seq -i $tf1 -j $tf2 -p blastn -D 1 -S 1 |");
    my @best;
    # save the longest alignment;
    while (<$run>) {
      next if /^\#/;
      my @row = split;
      next if $row[2] < $Min_PID; # require a minimum percent ID
      if ( ! defined $best[0] || $best[3] < $row[3] ) {
	@best = @row;
      }
    }
    if ( @best ) {
      $length_TIR = $best[3];
      $l_pos = $best[6]-1; # strings start at 0 not 1
      $r_pos = $best[8]-1;  # strings start at 0 not 1
      $longest_com_sub = substr($leftend,$l_pos,$length_TIR);
    }
  } else {
    for(my $start=0; $start <= length($leftend)-3;$start++) {
      for(my $len = 3; $len <=length($leftend)-$start; $len++) {
	my $sub = substr($leftend,$start,$len);
	my $pos = 1 + index($rightend,$sub);
	if($pos>=1 && length($sub) >= $length_TIR) {
	  $longest_com_sub = $sub;
	  $length_TIR = length($sub);
	}
      }
    }
    $l_pos = index($leftend,$longest_com_sub);
    $r_pos = index($rightend,$longest_com_sub);
  }
  next if ! $length_TIR;
  printf"TIR length: %d	'%s'\n",$length_TIR, $longest_com_sub;

  #make sure the nts before are expected TSD seqs
  printf "left: %d 		right: %d\n", $l_pos, $r_pos;
  my $leftTSD = substr($leftend,$l_pos-$TSD_length,$TSD_length);
  my $rightTSD = substr($rightend,$r_pos-$TSD_length, $TSD_length);

  printf "leftTSD: %s 		rightTSD: %s\n",$leftTSD,$rightTSD;

  if ( $length_TIR > $TIR_min_length &&
       $leftTSD eq $TSD_f && 
       $rightTSD eq $TSD_r ) {
      my $TE_start = $left_s + $l_pos-$TSD_length;
      my $TE_end   = $right_e - $r_pos + $TSD_length;
      my $full_TE = $dbh->seq($qid,$TE_start => $TE_end);
      printf "full: %s\n", $full_TE;
      printf "TSD: %s <--> %s\n", substr($full_TE, 0, $TSD_length), substr($full_TE, -1 * $TSD_length);
      my $TIR_right = (substr($full_TE, length($full_TE) - ($length_TIR + $TSD_length),
			      $length_TIR));
      my $TIR_left = substr($full_TE, $TSD_length, $length_TIR);
      printf "TIR:\n %s [fwd TIR %d..%d]\n %s [rev TIR revcom %d..%d]\n %s [rev TIR fwd]\n\n", 
      $TIR_left,
      $TSD_length, $TSD_length+ $length_TIR,
      revcom_as_string($TIR_right), length($full_TE) - ($length_TIR + $TSD_length), 
      length($full_TE) - ($length_TIR + $TSD_length) + $length_TIR,
      $TIR_right;
      
      my $TIR_right_index = rindex($full_TE, substr($full_TE, length($full_TE) - ($length_TIR + $TSD_length),
						    $length_TIR));
      my $TE_only = substr($full_TE, $TSD_length + $length_TIR,
			   $TIR_right_index - ($length_TIR + $TSD_length),
			   );
      printf "TE: %s \n", $TE_only;

      my $cased_full_TE = $leftTSD . lc($TIR_left) . $TE_only . lc($TIR_right) . $rightTSD;
      my $out = Bio::SeqIO->new(-format => 'fasta',
				-file   => ">$outdir/$Prefix\_$ct.fa");
      $out->write_seq(Bio::Seq->new(-id => "$Prefix\_$ct",
				    -desc => sprintf("%s:%d..%d TIR_length=%d TIR=%s leftTSD=%s rightTSD=%s",
						     $qid,$TE_start, $TE_end,
						     $length_TIR,
						     $TIR_left,
						     $leftTSD,$rightTSD),
				    -seq => $cased_full_TE));
      $ct++;
  }
}




